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SUMMARY 

This paper presents the general accuracy law that rules the MonteCarlo, ray-tracing algorithms 
used commonly for the calculation of the radiative entities in the thermal analysis of spacecraft. These enti- 
ties involve transfer of radiative energy either from a single source to a target (e.g., the configuration fac- 
tors). or from several sources to a target (e.g., the absorbed heat fluxes). In fact, the former is just a 
particular case of the latter. The accuracy model is later applied to the calculation of some specific radiative 
entities. Furthermore, some issues related to the implementation of such a model in a software tool are dis- 
cussed. Although only the relative error is considered through the discussion, similar results can be derived 
for the absolute error. 


INTRODUCTION 

Monte Carlo methods are often used in combination with ray-tracing algorithms to perform the 
radiative analysis of spacecraft (ref. [1]). Using this approach (MCFT from now, on), the radiative cou- 
plings between the faces of a model, as well as the external heat loads applied on these faces, are calcu- 
lated. Normally, these radiative values are passed to a thermal solver in order to produce the temperature 
predictions for the spacecraft model. 

While MCRT-based tools present some interesting advantages with respect to other methods, they 
also show a major drawback, which is the often large computational effort required to produce the radiative 
values. The results of a MCRT simulation are taken as an estimation of the actual values of the radiative 
entities. Since these results are of random nature, the accuracy of the estimation depends on the number of 
rays fired in the simulation. In general, the thermal engineer performs a trade off between the accuracy of 
the results and the computational effort required to achieve them. This paper presents models which allow 
the automatic accuracy control of the results in an efficient way. 

It is important to point out that the simulation inaccuracies only account for a pan of the uncer- 
tainty in the results of the thermal analysis. Other sources of error are the validity of the modelling assump- 
tions and the uncertainty of the data used in these models (ref. [2]). The “engineering judgement" shall be 
used to decide which level of accuracy in the simulation is sensibly required, especially when compared to 
the uncertainties in the other areas. Once the simulation's accuracy requirements have been established, 
their efficient achievement can be guaranteed by the techniques presented in tftis paper. 
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basic accuracy model 


. , . , ®ethods are based on the repetition of a gi ven random process, which reproduces the phys- 

°f inte , reSI For eacfl trial. values of the random parameters that play a role in the 

«f »h h UD \ ° rm y sam P ,ed from cumulative distribution functions, and a score T. representative 
e physic phenomenon, is tallied. T is a random variable that follows an arbitrary distribution: 


T- (E(T).V(T)) 


(EQ1) 


p h r,h E(T) . iS f e fC eCtati ° n ° f Tand V(T> varia "“ of T - estimation r* of the radiative parameter 

R is then calculated by averaging the scores over a large number of trials: 
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inS 1 fi ^ St0rC f ? r thC 11121 2,1(3 N is thc !otaI numbcr of trials. In the MCRT simulations. each trial 

J? f .J? ng 2 ray from source. Because of the random narure of the process, every simulation pro- 

rji n ated VaJ x, e ’ “ d therefore ^ estimation R* is a random variable itself. By the central 

of the disTrihut^n fn N rc *f onably Iar 2 e > R * is normally distributed, regardless of the actual shape 

of the distnbuuon for the basic random variable "R * 


R*-N(E(R*),V(R*)> 


(EO 3) 


where: 


E (R*) *p R . = 
V(R*)«o 2 r . 


E (T) = R 
V (T) 

3 
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By definition, the relative error of an estimated value r* is: 


(EQ 4) 
(EO 5) 
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(EQ 6) 


,ariabl “' “ is poss ‘ b,e 10 stK>w tt,at OTor ° f “* 


48 



E-N(E(E), V(E)) 


(EQ7) 


with expectation and variance given by: 


E(E)«ji £ = 0 


V(E )mj £ 


V(R «) 

R 2 


V (T) 
NR 2 


(EO 8) 
<EQ 9) 


Given the normal law followed by the simulation error, the probability o of having a relative error smaller 
than £ is: 


a«Prob (e <e) = erf 


and replacing (eq 9) into (eq 10) the fundamental accuracy model is derived: 



(EO 10) 
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V (T) 
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(EQ 11) 


This expression is very important for our purposes, because it provides a relationship between the 
accuracy e to be achieved and the computational effort (number of rays N E ) required for it, in terms of the 
radiative value itself, the variance of the basic random process used in the MCRT simulation and the confi- 
dence level a. Furthermore, (eq 11) shows that the accuracy achieved is inversely proportional to the 
square root of the number of rays fired in the simulation (i. e., to halve the relative error the number of rays 
has to be multiplied by four). 

For a given confidence level a, if N e rays are fired to ensure a level of accuracy of £. the variance 
of the relative error depends only on e: 



(EO 12) ■ 


That is, the variance of the relative error is directly proportional to the square of the level of accuracy spec- 
ified. Once the accuracy requirements are fixed, the variance of the errors associated to the estimation of 
different values do not depend on the actual values or on the variance of the basic random processes. This 
is an interesting property which is used later in the paper to introduce an efficient way to calculate recipro- 
cal couplings. 
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accuracy control for radiative parameters 

Configuration Factors 


.. . . con ^2 urau ° n (°[ v > ew ) favors represent the fraction of diffuse energy uniformly emined hy a 

radiauve face and reaching directly (j. e., without suffering any intermediate reflection) a target face. The 
procedure followed in the MCRT approach to determine the configuration factors is based on the repetition 
of a process consisung of the following steps: 


calculate a random emission point on the emitting radiative face 
• calculate a random emission direction 


utlculate the intersection between the ray fired from the emission point and directed along the emission 

direction, and the target radiative face 


• tally the ray if the intersection is not void 

Finally, the estimated value f* is calculated as the average of all the trials performed: 


1 V 
N ' X'i 


(EQ 13) 


I- 1 


Thus, this algorithm tallies the random variable T, which in this case represents the intersection 
between a randomly emitted ray and the target radiative face. T is a discrete random variable, which onl v 
can take two possible values: 0 if the intersection is void and 1 if the ray strikes the target radiative face. 
Since this is a rather simple process, the distribution function of the random variable T can be determined 
by only knowing the value of the configuration factor F, as shown in Figure 1. 

The expectation and variance of T can be calculated by applying the expressions used for discrete 
random variables: 


E(T) = £vP(, k ) 

k « I 
s 

V < T ) = X <t k ~E(T)) 2 .P(t k ) 

k- 1 


(EQ 14} 


(EO 15) 


fnT r ^th S ^ nUmbCr ° f ? SCretC Va,UCS which the random variablc «n take. Applying these expressions 
to T, with s - 2 corresponding to the only two possible values 0 and 1 , we obtain: 


E(T) = F 
V (T) = F- ( 1 - F) 

where F is the actual value of the configuration factor. 


(EQ 16) 
(EQ 17) 
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i urn 


Introducing these values into the general accuracy model given by (eq 1 1), the accuracy control law for 
the configuration factors can be derived: 



From this expression it can be appreciated that, in agreement with the experience, bounding the accuracy 
for small configuration factors is much more computationally demanding than doing so for laige values. 

Due to the fact that the variance of the basic random process depends exclusively on the value of 
the configuration factor, once this value F is fixed it is possible to specify the number of rays to be fired in 
order to achieve the desired accuracy. Figure 2 shows the computational effort (i. e. number of rays fired 
from the radiative face) required to achieve three different accuracy levels (relative errors of 1009c, 109r 
and \ %) for the whole range of configuration factor values. * 
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Configuration Factor F 


FIGURE 2. 


Computational effort vs configuration factor value for several levels of accuracy (confidence 
level of 99%) 


Radiative exchange factors 

The radiative exchange factors (REF) are defined as the fraction of diffuse energy uniformly emir- 
ted by a radiative face and being finally absorbed by a target face. Multi-reflection paths are included in 
this definition. These values can be derived from the configuration factors, following Gebhan's method 
(ref. [5]). Alternatively, the MCRT approach offers some advantages, taking into account the non-uniform 
nature of the radiati ve transfer exchange between radiative faces, as well as allowing the inclusion of spec- 
ular behaviour. The procedure followed for each radiative face is in this case: 

• calculate a random emission point on the emitting radiative face 

• calculate a random emission direction 

propagate the ray fired from the emission point and directed along the emission direction through the 
model. For the propagation take into account the radiative behaviour of the surfaces 

• tally the fraction of the original ray’s energy which is finally absorbed by the target face (random varia- 
ble T) 

Following this algorithm, the estimation of the REF value G is calculated by the expression: 


g* 



!- 1 


(EG 19) 
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The main difference with the configuration factors case is that it is not possible to know a priori the 
nature of the probability distribution function of the random variable T. Indeed. T can take different values 
depending on the path followed by the rays, and the number of possible paths grows significantly with the 
number of radiative faces in the model. Nevertheless, regardless of this fact the powerful central limit the- 
orem is valid and (eq 1 1) still applies. Therefore, the accuracy control model for REF can be expressed a.v 


N 


c 



(EQ 20) 


Thus, in the REF case, the computational effort depends explicitly not only on the value of the REF but 
also on the variance of the basic random process. 


Direct Solar Flux 


The heat flux emitted by the Sun and being intercepted by a given target radiative face can be cal- 
culated as: 


D = SC • A 1 


(EQ 21) 


where SC is the solar constant and A^* is the visible cross section area of the face. 

The MCRT method can also be used to calculate an estimation of D. Once the solar constant and 
the solar aspect ratio are known, the problem is reduced to find the visible cross section area of the radia- 
tive face. Although for non-oeduded planar faces this is a simple operation, a MCRT procedure can be fol- 
lowed whenever the faces are curved or shading effects exist: 

• calculate a random emission point on the radiative face 

• fire a ray from the emission point towards the sun 

• find whether the emission point “sees” the Sun. A discrete random variable H, taking only two possible 
values (0 if the Sun is not visible and 1 if it is) shall be tallied for this purpose 

• tally the cos 9 value (where 9 is the angle between the Sun direction and the face's normal vector at the 
emission point), only if the emission point is not occluded by any other part of the model 

An estimation of A x is then calculated as: 


a± * * jjf-A- X (cos0),-h, 
i - 1 


(EQ 22) 


where A is the radiative node area, h t is the value of the random variable H and (cos 0)[ is the value of the 
random variable cos 0. The subscript 1 refers to the 1-th trial. 

Applying the general accuracy model to this specific case, the following equation is obtained: 
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(EG 23) 


SOME IMPLEMENTATION ISSUES 
Parameter pre-estimation 


The application of the accuracy models to the estimation of radiative values presents the paradox 
of requiring as input the value whose calculation is the goal of the simulati on. Furthermore, the variance of 
the basic random process used to calculate' the estimation is not generally known a priori. These apparent - 
drawbacks can "be satisfactorily overcome by pre-estimating both the radiative value R and the variance of 
the basic process V(T), so that the accuracy models can be applied. The pre -estimated values can be 
obtained after a first batch of M rays is fired in the simulation. Indeed: 

M , 

R «t=X^ (EQ24) 

1- I 

V(T) «S 2 t * • j (t,-t) 2 (EQ 25) 

I- 1 


The number of rays M to be fired in order to pre-estimate R and V(T) shall be determined as a 
compromise between having reasonably accurate pre-estimations and not spending too much computa- 
tional effort in this previous phase of the algorithm. In practice, it has been checked that even with sample 
sizes which are small when compared to N c , the accuracy control based on the accuracy model produces 
excellent results. 

For instance. Figure 3 shows the histogram of the relative error associated to the estimation of a 
particular configuration factor, with a reference value of 0.01832. To produce the histogram, 1000 different 
simulations were performed, each of them using a pre-estimation sample size of 1000 rays. The accuracy 
model was then applied to ensure an accuracy of 3% with a level of confidence of 99%. The application of 
(eq 11) shows that these requirements are achieved by firing approximately 398000 rays in each simula- 
tion. The tails of the histogram, filled in black, show that only 10 out of the 1000 simulations performed 
had an error beyond the specified limits. This is in agreement with the confidence level used for this partic- 
ular case. 
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FIGURE 3. Histogram of the relative error associated to the estimation of a configuration factor with 
reference value 0.01832 


A particular problem arising when introducing the pre-estimated values into (eq 11) is that tills 
expression becomes singular for t * 0. A null pre-estimated value can indicate either that: 

• the actual radiative value is indeed null, in which case its actual variance will also be null, or that 

• the actual radiative value is so small that the size of the sample it is not large enough to provide a non- 
null pre-estimated value. 

In the first cast, no additional rays need to be fired. In the second, a very large number of rays are likely to 
be needed. Practical considerations impose a limit to this number, which for very small couplings might 
become computationally prohibitive. Generally speaking, once this limit is imposed it will not be possible 
to guarantee the accuracy of the radiative values below a threshold value. 


Enforcement of the reciprocity law 

In general, and for efficiency reasons, the software tools that implement MCRT methods do not 
calculate the couplings individually. Indeed, the couplings from one face to all the other faces in the model 
are normally calculated in one pass. Due to this fact, the couplings* line sum adds up to 1. However, the 
reciprocity law between couplings is in principle not satisfied, because of the statistical inaccuracy associ- 
ated to the estimations. 

Often, the manipulation of couplings that satisfy the reciprocity law is preferred, especially 
because of the reduction in memory requirements. The enforcement of the reciprocity law brings about a 
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Une sum which does not add up to 1. This imbalance is sometimes taken as a measure of the error in the 
estimation of the couplings. The enforcement of the couplings’ reciprocity is normally performed after 
both couplings have been independently calculated. This is not very efficient, since computational time is 
wasted to guarantee the specified level of accuracy for both the direct and inverse couplings. Furthermore, 
a systematic error is introduced by the enforcement method. 

On the other hand, the general error model shows that, given a reciprocal pair of faces, one of the 
asociated couplings is privileged in the sense that identical accuracy levels can be achieved with a smaller 
number of rays. This fact suggests the idea of estimating the non-privileged coupling by simply applying 
the reciprocity law to die one calculated via MCRT from the privileged face. Of course, this operation shall 
ensure that both coupling s errors meet the specified accuracy requirements. 

To apply this approach, the privileged face shall be identified. Therefore, the question to he 
answered is: given a reciprocal pair of couplings R,j and R Jt such that: 


where 'F is the reciprocity factor, and assuming that the accuracy requirements are respectively set to e 
and Cjj, which face shall.be used to fire rays from? lJ 

To identify the privileged face, let’s assume that the MCRT simulation is performed by firing rays 
from face 1 . Therefore, the coupling R^ is directly estimated via the simulation, while the couplin'* R, is 
estimated by enforcing the reciprocity law: e y 


r‘» = r* 

u u 

(EQ 27) 

r*.. 

r-.. a — 

(EQ 28) 


The reciprocal values estimated in such a way follow normal distributions: 


R V R % 


R's — N(E(R'«),V ( R* u ) } 


(EO 29) 
(EQ 30) 


Furthermore, it can be shown that the variance of the relative error associated to the 
(eq 28) is: 


estimation of R^ via 


V(E^) = V(Ej) 


V R « 


(EQ 31) 


On the other hand, if the estimation of Rjj was calculated by the MCRT procedure, by firing from 
face j the number of rays needed to meet the accuracy target Ej,, the variance of the relative error would he: 
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(E0 32) 


V(E S ) 


2 l erf" 1 (a) 
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Since the accuracy requirements must be also achieved by the enforced estimation R'.,, the variance of its 
error must be equal to the variance of the error obtained in the MCRT procedure: 


V(EJ = V(E t ) 

' i' 


(EQ 33) 


From this equation, it can be proved that i is the privileged face if the following condition holds: 


VHV) S'r-VfTV) 


(EQ 34) 


This expression relates the variance of the basic random processes used to estimate the couplings. If (eq 
34) is applied to the configuration factors case, it can be seen that the privileged face is the one with larger 
associated coupling. In the REF case, the privileged face cannot be determined by simply looking at the 
relative size of the reciprocal couplings, and (eq 34) shall be used instead. 

To clarify the interest of firing rays from the privileged face, let's present an example. For the 
model in Figure 4, reference values are available in the literature (ref. {6]). In particular, F 12 = 0.29176. 
Assuming we are interested in the calculation of the view factors with an accuracy of 5 9c, the most effi- 
cient way to proceed is to fire rays from face 1, which is the privileged one. This can be seen by checking 
the condition given by (eq 34). If the accuracy law for view factors (eq 18) is applied, it can be seen that 
approximately 6500 rays are needed to guarantee the accuracy requirements for both reciprocal couplings. 

If the same accuracy level had to be achieved by firing rays from face 2, roughlv 25000 rays would be 
needed. 



FIGURE 4. Model consisting of two perpendicular rectangular plates, with areas of 1 and 3 units 
respectively. The reciprocity factor between the configuration factors is in this case =3. 
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ACCURACY CONTROL FOR ABSORBED HEAT FLUXES 


The estimation of the absorbed heat flux on a radiative face differs of that of the simple radiative 
values presented above in that the absorbed value can be seen as the contribution of several source terms. 
Indeed, the heat flux striking the spacecraft follows different multi-reflection paths until it is finally 
absorbed by the face. Obviously, some multi-reflection paths will contribute more than others to the final 
absorbed flux. The purpose of this section is to apply the general accuracy model to this problem, taking 
into account the relevance of the mentioned radiative paths into the final value. The results that follow arc 
applicable to either solar or planetary (infrared and albedo) absorbed heat fluxes. 

A proper accuracy control for the absorbed heat flux is especially interesting when trying to quan- 
tify the absorbed heat flux on a spacecraft radiator. Because of the radiator's heat rejection requirements, 
the direct heat loads are generally small, and most of the absorbed flux reaches the radiator after several 
reflections. 

As previously stated, the heat flux <J> absorbed by a given radiative face can be expressed as the 
addition of a number of contributing terms 


♦-Z*i (EO 35) 

j- I 


where each <t>j is the flux being absorbed by the target face, with origin in the reflection of direct heat flux 
in the source face j. The summation is therefore extended to the n faces in the model which have non-null 
direct heat flux. 

The concept of source faces is clarified with the help of Figure 5. Assuming an idealised model 
consisting of 4 surfaces illuminated by the Sun, only faces 1 and 3 are source faces, as far as the calculation 
of the heat flux absorbed by face 4 is concerned. Face 2 does not see the Sun, and therefore it will not con- 
tribute with a source term in (eq 35). 
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FIGURE 5. 


face ?j 1 and 1 . 3 contribute to the heat flux absorbed by face 4. In this case, 
contributes through multi-reflections. 


iace i oniv 


au .i" . the a K CtUal ^ CRT simulation ' <t>j term is estimated by uniformly firing N. rays from face j 
by the tar^facr ran<J ° m vanable Ti ’ u e - amount of the energy carried by the ray that is finally absorbed 



I y 1 

n" * 2 


j i- 1 


The estimated value for the absorbed heat flux is therefore: 


(EQ 36) 


4 >* 



(EQ 37) 


and follows a normal law: 


<&*-N(E(<t>'),V (<!>*)) 


The variance of the relative error associated to the estimation is: 


(EQ 38) 
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(EO 39) 


V( E ). £v, 


r v £0 

N; 


This expression shows that the variance of the relative error, which is directly related to the accuracy- 
achieved by the simulation, depends on several terms, one for each source face. The question here is, given 
a fixed amount N of rays to be fired: 


N - i>, 

j - 1 


(EQ 40) 


how shall these N rays be distributed among the source faces so that (eq 39) is minimised? The solution to 
this problem can be found by applying standard techniques of non-linear programming. It can be shown 
that the optimal solution is given by: ~ 



•FTf~) 

0 

2yv ( T„) 

k» 1 



(EO 41) 


For this optimal distribution of the rays, the variance of the error becomes: 


V (E) 



(EQ 42) 


Replacing (eq ‘42) in the general error model, the accuracy control law for the estimation of the absorbed 
heat fluxes can be obtained: 



(EO 43) 


This multi-source model is the most general expression of the accuracy law for MCRT calcula- 
tions. In fact, the accuracy model previously obtained for the simple radiative entities is just a particular 
case of (eq 43), with n = 1. Indeed, one can regard the calculations of couplings as a mono-source phe- 
nomenon, being the unique source the radiative face emitting the radiation. Similarly, the Sun (or the 
planet) are the unique sources emitting radiation when the direct fluxes are calculated. 

The considerations about the pre-estimation of the values are also relevant for the implementation 
of this accuracy model. 
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CONCLUSIONS 


This paper has presented a model which can be used to control the accuracy of the radiative valuer 
estimated by using MCRT methods. The generality of the model allows its application to the radiative val- 
ues relevant to the thermal analysis of spacecraft. 

A prototyping activity has confirmed the validity of the assumptions used in the derivation of the 
accurac) models, as well as the feasibility of their implementation. The results obtained are in excellent 

agreement with the values predicted by the statistical models, even if pre-estimation of some parameters is 
required. 

^ scal * implementation of the error models is in progress, with the intention to enhance 
ESARAD, the radiative analysis software developed for the European Space Agency (ref. [7]). Significant 

improvements in terms of accuracy control, efficiency and performance are expected in relation to other 
radiauve codes currently available in Europe. 


NOMENCLATURE 

The following conventions have been followed in naming the random variables: 

Star indicates estimation of a radiative entity For example, F* represents the random variable "estima- 
tion of the view factor F\ 

Upper case is reserved to denote the distribution of the random variable, while lower case denotes a 
sample value from this distribution. For example, f* is the estimation of the view factor F, as provided 
by a single simulation run. 

Furthermore, the notation: 

R - N (E(R), V (R) ) 

mdicates that the random variable R follows a normal distribution with expectation E(R) and variance 
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